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ABSTRACT 

This  thesis  investigates  the  feasibility  of  micro- 
computer based  simulation  of  scalar  wave  propagation  in  var- 
ious media.  Models  for  lossless  media  and  media  with  a  loss 
coefficient  which  is  linear  in  frequency  have  been  coded  in 
FORTRAN  and  simulated  successfully  on  a  commercially  avail- 
able micro-computer,  with  simulation  times  less  than  thirty 
minutes.  The  spatial  impulse  responses  for  classical 
problems  using  square  and  circular-piston  excitation  are 
presented  graphically,  along  with  new,  innovative,  spatial 
excitation  source  shapes. 
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I.  INTRODUCTION 

Due  to  diffraction  effects,  acoustic  imaging  and  tissue 
characterization  techniques  require  information  about  the 
insonifying  wave  at  the  object's  location.  The  propagation 
of  a  pulsed  ultrasound  wave  with  arbitrary  temporal  and 
spatial  shape  is  not  well  understood,  and  requires 
development  of  computer-aided  analysis  and  simulation 
techniques. 

A  computationally  efficient  method  to  model  ultrasound 
propagation  from  a  planar  source  within,  a  medium  has  been 
developed  by  Guyomar  and  Powers  [Refs.  -1,2].  The  method  is 
applicable  to  lossless  media,  media  with  a  linear  frequency 
loss  coefficient  such  as  biological  tissue,  and  media  with  a 
quadratic  frequency  loss  coefficient  such  as  liquids  and 
gases.  Relying  on  Fast  Fourier  Transforms  (FFTs)  instead  of 
integral  solutions  that  require  geometric  interpretation 
makes  the  method  very  suitable  for  computer  implementation. 

Models  for  all  three  cases  have  been  previously 
developed,  coded  in  FORTRAN,  and  simulated  on  an  IBM  3  03  3 
mainframe  computer  [Ref.  3].  The  purpose  of  this  thesis  was 
to  investigate  the  feasibility  of  generating  micro-computer 
solutions  of  the  models  within  a  realistic  amount  of  time. 
As  the  average  execution  time  was  thirty  seconds  on  the 
mainframe  for  a  64  x  64  array  size,  a  time  limit  of  thirty 


minutes  was  set  as  a  realistic  goal.  While  the  difference 
between  thirty  seconds  and  thirty  minutes  may  appear 
unacceptable,  rapid  advances  in  micro-computer  technology  is 
narrowing  this  margin  appreciably. 

Since  the  previously  written  FORTRAN  programs  were  un- 
available, the  project  required  starting  from  the  beginning. 
A  first  attempt  was  made  using  a  commercially  available  pro- 
gram called  MATLAB.  MATLAB  is  an  array-oriented  program 
containing  built-in  functions  such  as  Bessel  and  Fast 
Fourier  Transforms,  as  well  as  the  ability  to  program  new 
functions.  Limitations  later  discovered  in  MATLAB  required 
starting  over,  programming  the  entire  project  in  Microsoft 
Fortran  Version  3.31.  The  result  is  a  stand-alone, 
user-interactive  program  allowing  customized  computer  runs 
tailored  to  a  specific  set  of  input  variables.  Only  the 
lossless  and  linear  loss  coefficient  cases  have  been  imple- 
mented, the  third  case  with  quadratic  loss  coefficient  is 
left  for  future  expansion.  The  goal  of  thirty  minutes  was 
obtained  by  code  optimization  and  exploiting  array  symmetry 
where  it  existed. 

The  program  was  tested  and  results  verified  using  clas- 
sical problems  such  as  a  square  or  circular-shaped  piston 
spatial  excitation.  Once  verified,  the  program  was  used  to 
compute  the  spatial  impulse  response  of  new  excitation 
shapes.  Examples  of  these  include  pyramid-shaped  pulsed 
sources  and  pulsed  linear  array  elements. 
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II.  PROBLEM  DESCRIPTION 

The  problem  solved  by  Guyomar  and  Powers  [Refs.  1,2] 
uses  the  geometry  of  Figure  1. 


OBSERVATION 
PLANE 


SOURCE 
PLANE 

Figure  1.  Source  and  Observation  Planes 


Given  a  separable  source  of  excitation 


v(x,y,0,t)  =  s(x,  y)T(t), 


(1) 


the  problem  is  to  find  the  acoustic  potential  <£(x,y,z,t)  at 
an  observation  point  located  a  distance  z  above  the  source 
plane.     A   rigidly   baffled   source   plane   and   linear, 


homogeneous  media  are  assumed.   It  has  been  shown  [Ref.  2] 
that  the  potential  is  given  by 


Hx>  y> 2> 0  =  5(x>  y)T(0 1 1 t 9fa  y,  *,  t) 


(2) 


The  *  indicates  convolution  over  the  variable  appearing 
directly  below  it  and  g(x,y,z,t)  is  the  impulse  response  or 
Green's  function  that  solves  the  wave  equation  and  meets  the 
applicable  boundary  conditions,  as  shown  in  the  block 
diagram  of  Figure  2. 


6(x.y.0.t) 


Propagation 

+ 

boundary   conditions 


g(x.y.z.t) 


^> 


Figure  2.  Block  diagram  of  impulse  response. 

In  the  spatial  domain,  0(x,y,z,t)  is  also  given  [Refs. 
4-8]  as 


<f>(x,  y,  z,  t)  =  T(t)  *t  h(x,  y,  z,  t) 


(3) 


where  h(x,y,z,t)  is  the  "spatial  impulse  response"  and  is 
equal  to 


h(x,  y,  z,  t)  =  s(x,  y)*x*  g(x,  y,  z,  t) 


(4) 


where  g(x,y,z,t)   is  the  Green's  function  (or  impulse  re- 
sponse) .   In  a  linear  system,  the  relationship  between  the 
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spatial  impulse  response,  h(x,y,z,t),  and  the  Green's 
function,  g(x,y,z,t),  is  illustrated  by  the  block  diagram  of 
Figure  3.   Using  the  two-dimensional   spatial  Fourier  trans- 


s(x.y)6(t) 


-O 


Propagation 

+ 

boundary   conditions 


htx.y.z.t)  = 
s(x.y)-**-g(x.y.z.t} 


Figure  3.  Block  diagram' of  spatial  impulse  response 
form  of  Equation  4,  tf>(x,y,z,t)  can  be  written  as 

4>(x,y,z,t)  =  T(t);T-l{sg} 


(5) 


where  the  tilde  notation  indicates  the  transform  of  the  spa- 
tial function. 

In  the  spatial  frequency  domain,  the  transform  of  the 
spatial  impulse  response  is 


h{x,y,z,t)  =  7      {sg} 


(6) 


The  convolutions  in  Equations  2  and  4  are  difficult  to 
implement  on  a  computer.  The  technique  presented  by 
Equations  3,  5,  and  6,  and  shown  graphically  in  the  block 
diagram  of  Figure  4,  reduce  two  of  the  convolutions  to 
multiplication. 


11 


From  Equation  6  we  see  that  if  we  interpret  s  as  the 
angular  spectrum  of  s,  then  tf  can  be  thought  to  be  the 
"propagation  transfer  function"  associated  with  propagation 


s(x.y)T(t) 


Propagation 

+ 

boundary   conditions 


-t> 


$  (x,y.z,t) 

=  T(t)  -*  h(x.y.z.t) 

t 

=  s(x.y)T(t)  *-*-*  g(x.y.z.t) 

xyt 

Figure  4.  Block  diagram  of  general  solution. 

in  the  medium.  This  propagation  transfer  function  behaves  as 
a  time-varying  spatial  filter  that  increasingly  attenuates 
the  higher  spatial  frequencies  of  the  source  as  time 
increases. 

Three  propagation  models  have  been  developed,  one  for 
each  type  of  media.  The  general  technique  followed  for  each 
model  begins  with  the  representative  wave  equation  for  the 
specific  medium  and  finding  the  Green's  function  (or  the  two 
dimensional  spatial  transform  of  the  Green's  function)  which 
solves  the  wave  equation.  Usinq  this  Green's  function  or 
its  spatial  transform,  the  spatial  impulse  response, 
h(x,y,z,t),  is  found  from  Equation  6.  From  this,  the  acous- 
tic potential  <f>  (x,  y ,  z  ,  t)  can  be  calculated  for  an  arbitrary 
plane  in  the  positive-z  half  space  using  Equation  5. 
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A.   CASE  1:   LOSSLESS  MEDIA 

The  results  of  the  transfer  function  approach  for  loss- 
less media  follow,  as  published  by  Guyomar  and  Powers  [Ref. 

1].  • 

The  wave  equation  is  the  Helmholtz  equation 


o         1  d26 
v2*-^=°-  (7) 


:2  dt2 
The  Green's    function   is 


8(ct  -  R)  .  . 

g(x,y,z,t)= — -— — ,  (8J 

Z7T  ii 


where 


R  =  sjx2  +  y2  +  z2  .  (9) 

For  this  problem  the  spatial  impulse  response  is 

h(x,  y,  z,  t)  =  f-— (10) 

ZTVrC 

Taking  the  spatial  transform  yields 

gpl   =  IsJ0  (py/cH*  -   22)  H(ct  -  z)  ,  (11) 

where   the   tilde   indicates   the   two-dimensional   spatial 
transform  and 

p  =  \Jfl  +  fi  ■  (12) 
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The  multiplicative  constants  can  be  dropped  since  the  re- 
sults will  be  normalized  to  maximum  values.  The  transfer 
function  for  propagation  in  lossless  media  is  then 


gpl  =  Jo  (pVcH2  -  z2)  H{ct  -  z) ,  (13) 


where  JQ  is  the  Bessel  function  of  the  first  kind  for  order 
zero.   For  programming  considerations  the  term 


z'  =  Vc2t2  -  z2  (14) 

is  calculated  as  a  unit  and  identified  as  ZPRIME  in  both 
this  paper  and  the  source  code. 

For  a  given  value  of  z,  one  can  calculate  the  value  of 
gpl  for  each  value  of  time,  then  take  the  inverse  spatial 
transform  of  the  product.  The  result  is  the  spatial  impulse 
response  of  Equation  8. 

B.   CASE  2:   LOSS  COEFFICIENT  LINEAR  IN  FREQUENCY 

For  media  which  exhibits  a  loss  coefficient  that  in- 
creases linearly  with  frequency,  the  results  of  Guyomar  and 
Powers  [Refs.  1,2]  follow. 

The  wave  equation  for  media  with  a  linear  loss  coeffi- 
cient is  the  telegrapher's  equation 

,    1  d24>         d<f>  ,  . 
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For  this  case  it  is  the  spatial  transform  of  the 
Green's  function  that  is  required.  It  has  been  found  by  an 
alternate  method  [Ref.  9]  to  be 

^    "jo  [Vp*-(A»cV4)v^^]  e-Ac^2H(ct  -  z)     for  p  <  Ac/2  ^ 

Jo  [jp2-{A2c2/A)y/cH2-z2\  e-Ac2t/2H(ct  -  *)     for  p  >  Ac/2 , 

where  IQ  is  the  Modified  Bessel  function  of  the  first  kind 
of  order  zero.  By  definition,  the  propagation  transfer 
function  is  equal  to  the  transform  of  the  Green's  function, 
so 


9p2{fx,fyjZ,t)  - 

I0  [vV  -  (A2c2/4)  Vc2t2  -  z2]  e-Ac2tf2H{ct  -  z)     for  p  <  Ac/2 


Jc 


y/p  -(A2c2/4)\/c2«2  -z2]  e-Ac2tl2H (ct  -  z)     for  p  >  Ac/2  . 


(17) 


The  transform  of  the  spatial  impulse  response  is 
Hfz,fv,z,t)  = 


s(/x,/»Ko    \/p2  -  (A2c2/4)  VcH2  -  z2    e-Actl2H{ct-  z)     for  p  <  Ac/2 


HfxJy)Jo    vV  -  (A2c2/4)  \/c2t2  -  z2    e~Actl2H(ct  -  z)     for  p  >  Ac/2. 


(18) 


The  algorithm  for  finding  the  spatial  impulse  response 
is  similar  to  the  lossless  case.  The  transform  is  found  of 
the  driving  function  s(x,y)  and  multiplied  by  gD2  evaluated 
at  each  spatial  frequency.  Taking  the  inverse  transform  of 
this  product  yields  the  required  spatial  impulse  response. 
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III.  PROGRAM  DESCRIPTION 

A.   PROGRAMMED  IN  FORTRAN 

Fortran  was  the  language  chosen  for  the  program  for 
three  reasons:  1)  the  availability  of  International 
Mathematical  and  Statistical  Library  (IMSL)  [Ref.  10]  rou- 
tines, specifically  those  used  to  calculate  Bessel  functions 
and  one-  or  two-dimensional  Fast  Fourier  Transforms,  2)  the 
portability  of  the  source  code  to  other  brands  of  micro- 
computers and  compilers,  and  3)  familiarity  of  the  source 
code  for  other  students  and  professionals.  Microsoft 
Fortran  Version  3.31  was  the  specific  brand  of  compiler 
used. 

B.   ADVANTAGE  OF  SYMMETRY 

Most  of  the  arrays  generated  by  the  program  are  sym- 
metrical. This  was  used  as  an  advantage  in  reducing  the 
number  of  calculations  required  and  therefore  the  execution 
time  of  the  program.  As  operation  of  the  program  is  de- 
scribed, the  areas  where  symmetry  has  been  used  are  identi- 
fied. Figure  5  is  provided  to  aid  in  these  explanations  by 
providing  a  "map"  of  an  array  base.  Figure  5  shows  a  64  x 
64  array  base  situated  on  the  X,Y  plane  and  divided  into 
four  quadrants.    Note  that  the  quadrants  are  unequal  in 
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size,   i.e.,   quadrant  II  includes  rows  1  through  3  3  and 
columns  1  through  3  3  for  a  33  x  33   "sub-array"   while  quad- 


1 


33 


64 


33 


64 


QUADRANT 
I  I 


QUADRANT 


I 


QUADRANT 
I  I  I 


OUADRANT 
IV 


>Y 


Y 
X 

Figure  5.  Base  array  configuration 


rant  I  includes  rows  1  through  33  and  columns  3  3  through  64 
for  a  32  x  32  "sub-array".  Similarly,  when  the  left  side 
(quadrants  II  and  III)  is  to  be  calculated  as  a  whole,  the 
calculations  cover  all  64  rows  and  columns  1  through  33. 
The  different  sizes  are  required  to  ensure  that  only  column 
33  contains  "peak"  information  and  also  that  element  (33,33) 
is  the  physical  center  of  the  array.   The  importance  of  Tihis 
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and  the  reason  why  it  is  required  will  be  explained  later. 
In-line  code  provides  the  required  "flipping"  actions  to 
expand  one  or  two  quadrants  into  a  full  64  x  64  array  of  in- 
formation. 

C.   BLOCK  DIAGRAM 

Figure  6  is  a  block  diagram  illustrating  the  infor- 
mation flow  through  the  program.  The  program  is  explained 
block  by  block  in  the  following  text. 

1.   User  input 

The  program  starts  by  obtaining  information  specific 
for  the  problem  to  be  run.  Using  on-screen  prompts,  the 
user  selects  a  problem  name,  filter  type,  excitation  shape 
and  size,  and  values  for  Z,  RHO,  and  MAXTIME.  If  filter  gp2 
is  selected  for  lossy  media,  the  user  is  further  prompted 
for  a  value  of  the  loss  coefficient  A. 

The  problem  name  assigned  must  be  eight  characters 
or  less  and  conform  to  MS-DOS  filename  specifications.  Once 
entered,  the  program  concatenates  three  different  extensions 
to  the  filename.  This  creates  the  names  of  the  three  output 
files  generated  by  the  program. 

The  first  file  created  is  <f ilename>. IMP,  where 
<filename>  is  the  filename  entered  by  the  user.  The  program 
stores  the  64  x  64  spatial  excitation  array,  s(x,y),  in  the 
default  directory  under  this  name  prior  to  taking  its  two 
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dimensional  FFT.    It  is  provided  should  the  user  want  to 
graphically  show  the  spatial  excitation  shape. 

The  second  file  created  is  <f ilename>.TXT.  This  is 
a  printable  text  file  which  summarizes  the  user's  inputs  for 
a  particular  problem. 


START 


USER 
INPUT 


Zpr  1  me 

VECTOR 


RH02 
ARRAY 


INPUT 
PULSE 


FILTER  g 

LOSSLESS 
MED  I  A 


FILTER 

LOSSY 
MED  I  A 


QP2 
LINEAR 


Figure  6.  Program  block  diagram. 

The  third  and  last  file  created  is  <f  ilename> .  DAT 
This  file  contains  the  final  64  x  64  array  which  is  the  spa- 
tial impulse  response  for  the  given  problem. 
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All  three  files  are  in  ASCII  text  format.  They  can 
be  printed  using  the  DOS  PRINT  command,  put  into  a  word 
processing  program  for  printing,  or  put  into  a  graphics 
program  for  conversion  into  a  picture.  All  figures 
generated  for  this  thesis  were  produced  by  translating  the 
<filename>. IMP  or  <f ilename>. DAT  file  to  MATLAB  format  using 
the  translation  program  supplied  with  MATLAB.  The  converted 
files  were  then  read  into  MATLAB  and  plotted  using  the  MESH 
command.  As  new,  improved,  graphic  routines  for  micro- 
computers are  introduced,  a  future  improvement  would  be 
incorporation  of  graphic  routines  into  the  program. 

As  the  program  is  written,  fourteen  different  source 
excitation  configurations  are  allowed,  nine  of  which  have 
been  implemented.  The  choice  of  fourteen  different  configu- 
rations was  arbitrary,  providing  a  good  selection  of  choices 
and  ease  of  future  expansion.  Each  of  the  remaining  five 
configurations  can  be  implemented  by  adding  its  name  to  the 
screen  formatting  code  and  writing  the  FORTRAN  code  to  gen- 
erate the  excitation  shape.  The  desired  source  excitation 
shape  is  selected  by  entering  the  selection  number  found 
next  to  the  shape's  description. 

The  variable  RHO  represents  the  radial  distance  in 
the  spatial  frequency  domain  as  shown  in  Equation  12 .  It  is 
identified  by  the  Greek  character  p  in  Equations  11  through 
13  and  16  through  18  but  will  be  referred  to  as  RHO  to  match 
the  source  code.    RHO  is  used  to  create  a  1  x  64  vector 
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called  RH01  which  in  turn  is  used  to  generate  a  64  x  64  ar- 
ray called  RH02.  The  vector  RH01  is  unique  since  it  must 
start  approximately  at  -RHO  and  increment  until  it  is  equal 
to  0.0  in  RH01(33).  From  here  it  continues  to  increase  un- 
til it  equals  +RH0  in  RH01(64).  Since  the  vector  is  of  even 
length,  it  cannot  be  formed  simply  by  dividing  the  entire 
range  of  RHO  (-RHO  to  +RHO)  by  64  and  using  this  value  as 
the  increment.  To  do  so  would  create  equal  values  in 
RH01(32)  and  RH01(33)  which  defeats  the  requirement  that 
RH01(33)  contain  singular  information,  in  this  case  the 
value  0.0. 

The  largest  value  of  time,  in  seconds,  for  the  time 
dimension  is  represented  by  the  variable  MAXTIME.  For  each 
problem,  time  starts  at  the  instant  that  the  wave  front 
first  arrives  at  the  observation  plane.  This  is  when  time  = 
Z/C,  where  C  is  the  sound  velocity  (1500  meters/second) ,  and 
increments  linearly  until  it  equals  MAXTIME.  To  prevent  a 
run-time  error,  the  user  is  prompted  for  a  value  of  MAXTIME 
which  is  greater  than  Z/C  seconds. 

The  height  of  the  observation  plane  above  the  X,Y 
plane  is  represented  by  Z  (in  meters)  .  It  is  typically  a 
small  number,  several  hundredths  of  a  meter,  or  can  be  set 
to  zero.  A  value  of  0.1  meters  (10  cm)  was  used  for 
debugging  the  program  and  producing  the  results  found  in  the 
Numerical  Simulations  section. 
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2.   Generate  ZPRIME  vector 

At  the  completion  of  user  input,  the  program  starts 

to  generate  the  required  vectors  and  arrays.   The  first  of 

these,  the  TIME  vector,  is  initialized  and  then  used  to 

generate  the  ZPRIME  vector.   Both  vectors  are  1  x  64  in  size 

and  represent  a  series  of  time  steps  starting  at  time  =  Z/C 

and  increasing  to  MAXTIME. 

The  program  contains  a  variable  called  STEP.    Its 

default  value  of  4   (which  can  be  changed)   sets  rows  1 

through  4  of  the  final  RESULT  array  to  zero  after  all 

calculations   have   finished.     This   simulates   the   step 

function  H(ct-z) .    Since  any  information  contained  in  the 

first   four   rows   is   lost,   these   values   need   not   be 

calculated.   For  this  reason,  TIME  starts  at  TIME (STEP+1) , 

or  TIME (5),  with  a  value  of  Z/C  seconds,  and  increments 

linearly . until  it  reaches  MAXTIME  in  TIME(64).   The  smooth 

linear  progression  of  time  from  Z/C  to  MAXTIME  is  obtained 

using  the  following  code  segment: 

INC  =  (MAXTIME  -  TS) /REAL ( NROWS -STEP- 1) 
DO  XX  I  =  (STEP+1) , NROWS 

TIME (I)  =  TS  +  (I-STEP-1)  *  INC 
CONTINUE 

where  TS  =  Time  Start  =  Z/C,  NROWS  =  64,  and  STEP  =  4. 

As  an  example  let  MAXTIME  =  1.0E-4  seconds  and  Z  = 

0.1  meters  as  determined  by  the  user.    The  quantity  TS 

equals  6.666667E-5  seconds  and  this  value   is  stored  in 

TIME  (5).    Each  element  of  TIME  is   increased  by  INC  = 
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5.649718E-7  until  it  equals  MAXTIME  in  TIME(64).  The 
elapsed  time  for  this  particular  simulation  is  MAXTIME  -  TS 
=  3.333333E-5  seconds. 

The  ZPRIME  vector  is  created  using  values  found  in 
the  TIME  vector  as  shown  in  Equation  14.  Each  value  of  TIME 
is  squared  and  multiplied  by  C2 .  Subtracting  Z2  from  this 
product  and  taking  the  square  root  produces  a  ZPRIME  ele- 
ment. Continuing  with  the  previous  example,  two  sample 
ZPRIME  elements  are 

ZPRIME(6)  =  SQRT(TIME(6)2  *  C2  -  Z2)  =  1.304644E-2 

ZPRIME(7)  =  SQRT(TIME(7)2  *  C2  -  Z2)  =  1.848934E-2 
The  algorithm  loops  until  ZPRIME (STEP+1)  through  ZPRIME (64) 
have  been  calculated. 

3 .   Generation  of  RH02  array 

The  next  significant  step  is  formation  of  the  RH02 
array.  First  the  RH01  vector  is  created  by  calculating  an 
increment  and  base  value  using  the  code 

INC  =  REAL (RHO) /REAL (NCOLS/2-1) 

BASE  =  REAL ( -RHO 1)  -  INC 
Assuming  a  value  of  200  for  RHO,  INC  =  6.451613  and  BASE  = 
-206.451613.  This  results  in  a  RHOl  vector  starting  with 
-206.451613  in  RHOl(l),  with  each  successive  element  incre- 
mented by  INC.  This  method  ensures  a  value  of  0.0  in 
RH01(33),  -200.0  (-RHO)  in  RH01(2),  and  +200.0  (+RHO)  in 
RH01(64).  As  with  the  TIME  vector,  simply  setting  the  RHOl 
vector  to  range  from  -RHO  in  RHOl(l)  to  +RHO  in  RHOl  (64) 
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would  result  in  identical  values  in  RH01(32)  and  RH01(33). 
This  would  not  meet  the  requirement  that  RH01(33)  contain 
singular  information. 

After  completing  the  RH01  vector,  its  values  are 
used  to  generate  the  RH02  array.  First  the  RH01  vector 
values  are  placed  in  row  33  of  RH02  and  its  transpose  in 
column  33.  The  program  then  fills  in  the  blank  array 
elements  in  quadrants  II  and  III  using  the  Pythagorean 
theorem  with  the  values  found  in  the  respective  row  33  and 
column  3  3  positions. '  Continuing  with  the  example  using  RHO 
=  200,  the  value  -206.451613  would  be  found  in  elements 
RH02(3  3,1)  and  RH02(1,33)  (column  major  format).  Using  the 
Pythagorean  theorem,  element  RH02(1,1)  would  receive  the 
value  291.966671.  This  represents  the  value  of  RHO  as 
calculated  radially  from  the  central  position  of 
RH02(33,33). 

Since  the  RH02  array  is  symmetrical,  only  quadrants 
II  and  III  (left  side)  must  be  calculated.  Although  the 
time  saved  here  is  small,  the  algorithm  using  these  values 
later  in  the  program  needs  only  the  values  found  in  the  left 
half.  Therefore  the  right  side  values  are  not  calculated. 
4 .   Spatial  excitation  generation 

Generating  the  spatial  excitation  array  is  the  next 
step.  The  first  five  selections  represent  single  piston 
spatial  excitations  centered  about  position  (33,33)  on  a  64 
x  64  base  array  called  PULSE.   The  five  shapes  available  are 
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the  square  piston,  circle  piston,  gaussian,  pyramid,  and 
truncated-pyramid.  All  five  shapes  require  the  user  to 
specify  the  width  (or  diameter)  of  the  excitation.  The 
value  of  width  must  be  an  odd  integer  to  ensure  symmetry 
about  PULSE(33,33)  ,  and  less  than  or  equal  to  63  since  the 
base  array  size  is  64  x  64.  Small  values  between  9  and  15 
provide  excellent  results. 

The  square  and  circle  selections  create  a  simple 
spatial  excitation  with  amplitude  =  1.0.  A  resolution 
problem  exists  with  the  circular  excitation.  Simply,  it  is 
impossible  to  represent  a  circular  shape  using  a  small 
number  of  square  elements.  The  larger  the  diameter,  the 
smoother  the  outer  surface  becomes, "but  with  a  64  x  64  array 
size  the  range  of  available  diameters  is  limited.  The  only 
solution  is  an  increase  in  array  size  to  128  x  128  or 
larger.  This  is  an  excellent  area  for  future  expansion  as 
micro-computers  become  faster  and  less  memory  limited. 

The  Gaussian  selection  creates  a  circular-shaped 
excitation  with  amplitude  following  a  gaussian  distribution. 
The  spatial  excitation  is  truncated  at  the  user  specified 
diameter,  and  its  width  measured  at  the  1/e  point.  This 
excitation,  as  well  as  the  pyramid  and  truncated-pyramid, 
represent  non-piston  sources  with  spatial  variations  in 
amplitude. 

The  pyramid  creates  a  four-sided  tapered  spatial 
excitation  with  a  peak  amplitude  of  1.0  at  the  center  and 
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tapering  off  until  it  equals  0.0  at  the  base.  The  size  of 
the  base  is  determined  by  the  input  value  of  width.  For 
example,  a  pyramid  excitation  of  width  11  would  have  its 
center  at  PULSE (33 ,33)  and  occupy  the  region  bordered  by 
rows  29  through  37  and  columns  29  through  37.  The  amplitude 
would  decrease  linearly  and  equal  zero  at  the  outer  edges  of 
the  above  region.  The  truncated-pyramid  is  similar  except 
that  the  amplitude  tapers  off  at  a  slower  rate  until  it 
equals  0.0  at  the  array *s  outer  edges.  At  the  specified 
value  of  width,  the  pyramid  is  truncated,  resulting  in  a 
"square  house"  shape  with  tapered  roof  and  square  sides. 

The  remaining  selections  all  represent  linear  arrays 
of  five  elements  each.  An • individual  element  size  of  5  x  5 
located  on  10  unit  wide  center  points  was  selected  for 
convenience.  Shape  selections  for  the  individual 
excitations  are  square,  circular,  and  gaussian.  Two 
additional  arrays  also  use  5  square  elements,  but  the 
amplitude  of  the  elements  are  stepped  or  parabolic  in  shape. 
The  stepped  array  has  a  center  element  of  amplitude  =  1.0, 
two  adjacent  elements  of  amplitude  =  2/3,  and  two  outer 
elements  of  amplitude  =  1/3.  The  five  elements  thus  create 
a  stepped  appearance:  1/3,  2/3,  1.0,  2/3,  1/3.  The 
parabolic  array  is  similar  to  the  stepped  array,  except  the 
amplitudes  follow  a  parabolic  shape. 

After  the  spatial  excitation  is  generated,  the  two 
dimensional  Fast  Fourier  Transform  is  taken  using  the  IMSL 
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routines  FFT3D  and  FFTCC.  The  routine  FFT3D  computes  the 
FFT  of  a  complex-valued  two-dimensional  array  by  passing 
first  the  columns  and  then  the  rows  as  individual  vectors  to 
the  routine  FFTCC  which  computes  the  FFT  of  the  vector  argu- 
ment. The  array  returned  by  FFT3D  is  complex  and  has  its 
peak  information  located  at  the  outer  four  corners  of  the 
array,  with  the  largest  value  in  PULSE (1,1).  The  array  is 
shifted  to  bring  this  peak  information  to  the  center  using 
the  subroutine  SHIFT.  This  subroutine  swaps  quadrant  I  with 
quadrant  III  and  quadrant  II  with  quadrant  IV,  thus  trans- 
ferring the  largest  value  from  element  PULSE (1,1)  to  element 
PULSE(33, 33) .  This  is  why  it  is  so  important  to  arrange  all 
arrays  so  their  central  value  is  in  element  (33,33). 
5.   Filter  apl  for  lossless  media 

From  this  point  the  program  branches  based  on  the 
filter  type  selected.  The  algorithm  for  the  lossless  case, 
using  filter  gpi,  is  to  take  a  value  of  ZPRIME,  multiply  the 
RH02  array  by  this  value  and  pass  each  product  as  the 
argument  to  the  IMSL  subroutine  MMBSJN,  where  MMBSJN 
calculates  the  Bessel  function  of  the  first  kind  of  order 
zero  for  a  real  number  argument.  The  result  returned  from 
the  Bessel  function  is  multiplied  by  the  respective  element 
in  array  PULSE  and  stored  in  array  WORK.  The  next  value  of 
ZPRIME  is  selected  and  the  calculations  again  performed. 

As  this  portion  of  the  program  is  the  most  time 
consuming,  two  specific  actions  were  taken  to  optimize  the 
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code.  As  previously  explained,  the  number  of  rows 
identified  by  STEP  are  zeroed  in  the  final  array.  Because 
of  this,  calculations  start  at  row  (STEP+1) ,  thus  saving  the 
calculation  of  4  rows  for  each  iteration.  Also,  since  the 
arrays  RH02  and  PULSE  are  symmetrical  in  all  four  quadrants 
about  element  (33,33),  only  quadrant  II  needs  to  be 
calculated.  With  the  default  64  x  64  array  size  and  STEP  = 
4,  only  59  (64-STEP-l)  subarrays  of  size  33  x  33  must  be 
calculated  for  a  total  of  64,251  array  elements.  This 
compares  favorably  with  the  262,144  array  elements  that 
would  require  calculation  if  the  entire  64  x  64  array  was 
calculated  64  times.  It  is  important  to  note  that  the 
single  most  time-consuming  process  is  the  calculation  of  the 
Bessel  function. 

After  each  33  x  33  subarray  is  calculated,  inline 
code  is  used  to  fill  the  entire  WORK  array.  Quadrant  II  is 
flipped  to  fill  quadrant  III,  then  the  entire  left  side  is 
flipped  to  fill  the  right  side. 

At  this  point  the  inverse  two-dimensional  Fast 
Fourier  Transform  of  the  WORK  array  needs  to  be  taken.  But 
instead  of  passing  the  entire  array  as  an  argument  to  a  two 
dimensional  inverse  FFT  subroutine,  the  array  is  processed 
one  column  at  a  time.  Using  the  IMSL  subroutine  FFT2C, 
which  computes  the  inverse  FFT  (or  FFT)  of  a  complex-valued 
vector,  each  of  the  64  columns  are  passed  individually  as  a 
vector.    After  all  the  columns  are  transformed,  only  one 
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row,  row  33,  is  passed  as  an  argument.  Assuming  a  symmetri- 
cal transducer,  row  33  contains  information  from  the  central 
line  across  one  axis  of  the  observation  plane.  The  inverse 
transform  of  this  row  represents  the  field  values  along  this 
line  for  the  specific  value  of  ZPRIME,  and  it  is  this 
information  that  is  used  to  sequentially  build  the  final 
array  called  RESULT.  The  59  iterations  completed  by  the 
above  algorithm  will  each  supply  one  row,  corresponding  to 
an  instant  of  time,  in  array  RESULT.  Specifically,  with  the 
default  value  of  STEP,  rows  5  (STEP+1)  through  64  (NROWS) 
are  filled. 

After   completing   row   64,   the   number   of   rows 
identified  by  STEP  are  set  to  zero  as  they  do  not  contain 
any  information.   The  entire  RESULT  array  is  then  saved  to 
disk,  and  the  program  terminates  normally. 
6.   Filter  gp2  for  lossy  media 

In  the  second  case,  using  filter  gp2  for  media  with 
linear  loss  coefficient,  the  algorithm  has  two  additional 
steps.  The  program  implements  Equation  17  and  therefore 
must  test  each  argument  for  RHO  greater  than  AC/2.  For 
RHO  >  AC/2,  the  argument  is  passed  to  the  IMSL  Bessel 
function  subroutine  MMBSJN  as  in  the  lossless  case.  But,  if 
RHO  <  AC/2,  the  argument  is  instead  passed  to  the  IMSL 
subroutine  MMBSIN  which  calculates  the  modified  Bessel 
function  of  the  first  kind  of  order  zero.  The  result 
returned  from  either  Bessel  function  is  multiplied  by  the 
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respective  element  of  PULSE.  Then,  as  shown  in  Equation  17, 
the  product  is  multiplied  by  an  exponential  term  prior  to 
storing  in  array  WORK.  This  is  the  second  difference  in  the 
lossy  case. 

The  algorithm  loops,  as  in  the  lossless  case,  to 
fill  the  entire  RESULT  array.  After  saving  the  RESULT  array 
on  disk,  the  program  terminates  normally. 
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IV.  NUMERICAL  SIMULATION 

After  careful  debugging  and  testing  on  a  point  by  point 
basis,  the  program  was  used  to  calculate  the  spatial  impulse 
response  for  various  source  excitations.  For  each  excita- 
tion shape,  both  the  lossless  and  lossy  case  were  calculated 
and  the  results  plotted.  The  plots  show  the  amplitude  of 
the  wave  plotted  against  cross-direction  and  time  at  an  ob- 
servation plane  located  above  the  X,Y  plane.  In  all  calcu- 
lations the  height  of  the  observation  plane,  Z,  is  10  cm 
above  the  source  plane.  The  value  of  A  used  for  the  lossy 
cases  was  8.0E-3.  This  value  was  selected  as  a  compromise 
and  arrived  at  using  trial  and  error.  Smaller  values  of  A 
caused  only  minor  variations  in  shape  while  larger  values 
caused  gross  distortion.  All  problems,  unless  noted  other- 
wise, are  calculated  for  a  MAXTIME  of  1.5E-4  seconds.  Thus, 
with  a  Z  value  of  10  cm,  each  result  runs  from  a  starting 
time  =  Z/C  =  6.666667E-5  seconds  to  1.5E-4  seconds,  a  dura- 
tion of  8.333333E-5  seconds. 

As  previously  explained,  all  plots  were  obtained  using 
the  MATLAB  MESH  command.  Two  limitations  exist  with  the 
MESH  routine  that  reguire  further  explanation.  The  first 
limitation  is  the  way  MATLAB  sizes  the  plot.  MATLAB  uses 
input  array  values  to  fill  a  window  of  pre-determined  size. 
This  fixes  the  height  of  the  plot  so  an  excitation  with 
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amplitude  equal  to  1.0  will  plot  identically  to  an 
excitation  with  amplitude  other  than  1.0.  The  second 
limitation  builds  on  the  first:  MATLAB  does  not  show 
numerically  scaled  axes  when  plotting  an  array  of  data.  If 
it  did,  the  different  amplitudes  for  identically  shaped 
plots  would  be  apparent. 

An  alternative  way  to  obtain  numerical  information  was 
found.  Using  the  PLOT  command,  MATLAB  has  the  ability  to 
plot  a  vector  in  the  X,Y  domain  with  numerically  scaled 
axes.  Putting  an  array  of  data  into  PLOT  creates  a  side 
view,  plotting  successive . "slices"  of  the  array  with  time 
increasing  in  the  positive  X  direction.  Using  this  tech- 
nique, arplitude  values  could  be  obtained  for  a  given  plot. 
Since  shape  is  the  primary  interest,  amplitude  information 
obtained  in  this  manner  was  adequate.  Three  plotted  arrays 
(Figures  9,  11,  and  15)  are  included  for  illustration. 

Execution  times  were  all  less  than  the  stated  goal  of 
thirty  minutes,  with  an  average  of  twenty-five  minutes. 
These  times  were  obtained  using  an  AT&T  63  00  micro-computer 
operating  at  8  MHz.  An  8087  numerical  co-processor  chip  was 
installed  and  used  to  increase  performance.  Similar 
performance  can  be  obtained  with  fast  IBM  XT  compatibles  or 
AT  class  micro-computers.  The  next  generation  of 
micro-computers  utilizing  the  Intel  80386  microprocessor  are 
the  logical  choice  for  continued  program  development. 
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A.   CLASSICAL  EXAMPLES 

The  first  two  source  shapes  presented  are  the  square 
piston  and  circular  piston. 


Figure  7.  Square  piston  spatial  excitation. 

Figure  7  shows  graphically  the  source  excitation  shape 
for  a  square  transducer  of  width  11  and  amplitude  1.0.  The 
excitation  is  centered  over  the  base  array's  central  ele- 
ment, PULSE (33, 33) .  In  Figure  8  the  calculated  spatial  im- 
pulse response  for  the  square  transducer  in  lossless  media 
is  shown,  as  is  the  array's  orientation  with  time  and  space. 
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This  orientation  of  time  and  space  applies  to  all  spatial 
impulse  response  plots.  The  origin  of  the  time  axis  starts 
at  Z/C  where  the  potential  is  known  to  be  a  replica  of  the 
excitation  source  shape  for  lossless  media.  As  time 
progresses,  the  potential  is  reduced  until  it  forms  two  out- 
ward traveling  "tails".   As  expected,  the  "tails"  are  square 


time 


Figure  8.  Field  for  square  transducer  (lossless  case). 

to  match  the  spatial  excitation  shape,  and  slowly  attenuate 
with  time  in  the  lossless  case.  Numerical  information  can 
be  obtained  from  Figure  9  which  is  the  side  view  of  the 
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plot  in  Figure  8.  The  step  function  is  clearly  visible  as 
is  the  jump  in  potential  to  1.0  in  row  5.  From  this  point 
the  potential  peaks  at  a  value  slightly  greater  than  1.0  and 
then  drops  rapidly.  At  infinity  the  tails  should  reduce  to 
zero. 
1..  2. 
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Figure  9.  Side  view  for  lossless  case. 

In  Figure  10  the  same  excitation  has  been  analyzed  for 
the  lossy  case.  Comparing  the  shapes  in  Figures  8  and  10 
reveals  two  distinct  differences.  The  "tails"  are  no  longer 
square  but  have  taken  on  a  triangular  appearance.    Larger 
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values  of  A  will  produce  a  sharper  triangular  shape.  The 
second  major  difference  concerns  the  region  between  the 
tails.  In  the  lossless  case  this  region  was  at  zero  poten- 
tial. In  the  lossy  case,  the  region  does  not  approach  zero 
but  instead  fills  in,  developing  a  shape  of  its  own. 


Figure  10.  Field  for  square  transducer  (lossy  case). 

Another  major  difference  is  discovered  when  comparing 
Figures  9  and  11.  Where  the  initial  amplitude  for  the  loss- 
less case  was  1.0,  the  amplitude  of  the  spatial  excitation, 
the  initial  amplitude  for  the  lossy  case  has  been  reduced  to 
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approximately  0.30.    This  reduction  in  amplitude  was  ob- 
served for  all  test  cases  with  a  linear  loss  coefficient. 

The  results  obtained  using  a  circular  transducer  are 
similar  to  those  obtained  with  the  square  transducer.  The 
spatial  excitation  for  a  circular  transducer  of  width  15  and 
amplitude  1.0  is  shown  in  Figure  12.   Here  a  width  of  15  was 
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Figure  11.  Side  view  for  lossy  case. 

used  to  improve  the  resolution  of  the  circular  excitation. 
The  spatial  impulse  response  using  this  source  is  shown  in 
Figure  13.  Comparison  of  Figures  8  and  13  illustrates  the 
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differences  between  the  spatial  excitation  shapes.  The  ma- 
jor difference  is  the  tails,  which  now  have  a  circular 
cross-sectional  shape.  Some  minor  variations  can  also  be 
observed  in  the  large  initial  response  early  in  time. 
Figure  14  shows  the  spatial  impulse  response,  using  the 
circular  excitation  of  Figure  12,  for  the  lossy  case.  The 
previous   comments  for  the  square   transducer  concerning  the 


Figure  12.  Circular  piston  spatial  excitation. 

difference  in  amplitudes  and  shape  of  the  "tails"  for  the 
lossless   and   lossy   cases   also   apply   to   the   circular 
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transducer.  Figure  15  is  a  side  view  plot  of  Figure  14 
showing  the  reduction  in  amplitude  is  also  present  with  the 
circular  transducer. 

The  results  obtained  when  using  the  square-piston  and 
circular-piston  spatial  excitation  for  both  lossless  and 
lossy  media  compare  favorably  with  those  obtained  by  Guyomar 
and  Powers  [Refs.  1-3,9].  Based  on  these  results,  it  was 
determined  that  the  program  correctly  implemented  the 
propagation  models,  and  therefore  could  be  used  to  compute 
the  spatial  impulse  response  using  new  types  of  excitation. 


Figure  13.  Field  for  circular  transducer  (lossless  case) 
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B.   NEW  EXCITATION  SHAPES 

The  program  has  to  ability  to  analyze  any  excitation 
source  shape  that  will  fit  on  the  64  x  64  base  array.  Two 
new  single  excitation  shapes  and  four  multiple  element  exci- 
tation shapes  are  in  the  program.   The  results  for  three  of 


Figure  14.  Field  for  circular  transducer  (lossy  case). 

these,  the  pyramid,  the  5  element  linear  array  with  constant 
amplitude,  and  the  5  element  linear  array  with  stepped  am- 
plitude are  included. 
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A  pyramid  shaped  spatial  excitation  of  width  11  and  am- 
plitude 1.0  is  shown  in  Figure  16.  The  computed  spatial  im- 
pulse response  for  this  excitation  in  lossless  media  is 
shown   in  Figure  17.   Since  this  particular  waveshape  has  a 
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Figure  15.  Side  view  for  lossy  case. 

low  spatial  frequency  content,  the  shape  of  the  wave  stays 
much  the  same  for  small  values  of  time.  Figure  18  shows  the 
spatial  impulse  response  for  the  same  pyramid  excitation  in 
lossy  media. 
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Figure  19  shows  the  5-element  linear  array  configura- 
tion. For  this  array  the  elements  are  each  5x5  with  their 
centers  positioned  10  units  apart.  This  produces  a  space  of 
5  units  between  each  element.  The  calculated  spatial  im- 
pulse response  for  this  excitation  is  shown  in  Figure  20. 


Figure  16.  Pyramid  shaped  spatial  excitation. 

The  shape  for  each  individual  excitation  is  similar  to  the 
single  excitation  response  in  Figure  8  until  the  tails  start 
to  spread  out  and  interfere  with  each  other.   In  the  lossy 
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case,  Figure  21,  the  pattern  is  similar  but  with  increased 
interference  in  the  region  between  the  individual  responses. 


Figure  17.  Field  for  pyramid  excitation  (lossless  case). 

The  spatial  excitation  for  the  last  test  case  presented 
is  shown  in  Figure  22.  This  is  a  5-element  linear  array 
similar  to  Figure  19  but  with  stepped  amplitude.  The  center 
element  has  amplitude  of  1.0  as  before.  The  two  elements 
adjacent  to  the  center  have  amplitudes  of  2/3,  and  the  outer 
two  elements  have  amplitudes  of  1/3.  This  particular  case 
is  presented  to  show  the  flexibility  of  the  program.   Figure 
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2  3  shows  the  spatial  impulse  response  in  lossless  media. 
The  results  are  similar  to  those  of  Figure  20  except  for  the 
dominance  of  the  center  element.  Figure  24  shows  the  spa- 
tial impulse  response  for  the  same  excitation  in  lossy  me- 
dia. 


Figure  18.  Field  for  pyramid  excitation  (lossy  case) 
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Figure  19.  Excitation  for  five  element  linear  array. 
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Figure  20.  Field  for  five  element  linear  array 

(lossless  case) . 
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Figure  21.  Field  for  five  element  linear  array 

(lossy  case) . 
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Figure  22.  Excitation  for  five  element  linear  array 
with  stepped  amplitude. 
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Figure  23.  Field  of  five  element  linear  array  with 
stepped  amplitude  (lossless  case) . 
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Figure  24.  Field  for  five  element  linear  array  with 
stepped  amplitude  (lossy  case) . 
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V.  SUMMARY 

This  thesis  has  investigated  the  ability  to  perform 
complex  simulations  of  scalar  wave  propagation  through  loss- 
less media  and  media  with  a  loss  coefficient  linear  in  fre- 
quency on  a  micro-computer.  A  program  was  written  in 
Microsoft  Fortran  V3.31  which  allows  the  user  to  select  or 
specify  the  initial  conditions  and  media  type,  then  calcu- 
lates the  resulting  spatial  impulse  response.  The  set  goal 
of  thirty  minutes  for  each  simulation  was  achieved  through 
code  optimization  and  partial  array  calculation  permitted  by 
symmetry . 

After  initial  program  verification,  the  spatial  impulse 
responses  for  both  square-piston  and  circular-piston  spatial 
excitation  were  computed.  The  results  obtained  agreed  with 
previously  accepted  results.  The  spatial  impulse  response 
was  then  calculated  for  three  new  spatial  excitation 
sources,  a  pyramid-shaped  excitation,  a  five  element  linear 
array  excitation  with  constant  amplitude,  and  a  five  element 
linear  array  excitation  with  stepped  amplitude.  The 
computed  spatial  impulse  response  for  all  simulations  have 
been  plotted  and  are  presented  for  visual  interpretation. 

Future  development  should  concentrate  in  two  areas. 
First,  the  third  model  for  media  with  a  loss  coefficient 
which  is  quadratic  in  frequency  must  be  added  to  the 
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program.  This  will  provide  all  the  tools  necessary  for 
analysis  of  wave  propagation  through  air,  liquid,  gas,  and 
biological  tissue.  Second,  the  base  array  size  must  be  in- 
creased beyond  the  present  64  x  64.  An  array  size  of  64  x 
64  proved  adequate  for  program  development  and  testing  but 
increased  size,  and  therefore  increased  resolution,  is  nec- 
essary for  serious  analysis  and  evaluation.  The  limitation 
in  all  presently  available  micro-computers  which  restricts 
array  size  is  processing  speed.  For  example,  increasing  ar- 
ray size  to  128  x  128  would  increase  execution  time  by  a 
factor  of  four,  requiring  almost  two  hours  for  each  simula- 
tion run.  Transferring  program  development  to  one  of  the 
newer  micro-computers  utilizing  the  Intel  80386  micro- 
processor will  allow  larger  array  sizes  while  maintaining 
reasonable  execution  times. 
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APPENDIX:   FORTRAN  SOURCE  CODE 

This  appendix  contains  the  FORTRAN  source  code  for  the 
program  developed  in  the  thesis.  This  program  was  compiled 
and  run  using  Microsoft  Fortran  Version  3.31. 


$LARGE 

SNOFLOATCALLS 

C 

n  ****************************************************************** 

C  COMMENTS : 

C 

C  THIS  PROGRAM  WAS  WRITTEN  BY  LT .  TIMOTHY  MERRILL,  USN 

C 

C  IT  IS  WRITTEN  IN  MICROSOFT  FORTRAN  V3.31  AND  IS  CAPABLE  OF  RUNNING 

C  ON  ANY  IBM  OR  COMPATIBLE  MICRO-COMPUTER  WITH  LITTLE  OR  NO 

C  MODIFICATION.   THE  USE  OF  A  8087  NUMERICAL  CO-PROCESSOR  IS  HIGHLY 

C  RECOMMENDED  TO  REDUCE  EXECUTION  TIME. 

C 

C  TO  REDUCE  EXECUTION  TIME  THE  PROGRAM  IS  WRITTEN  TO  USE  A  64  X  64 

C  BASE  ARRAY  SIZE,  AND  IS  WRITTEN  AS  SINGLE  PRECISISON.   FUTURE 

C  MODIFICATION  TO  A  128  X  128  (OR  LARGER)  ARRAY  SIZE  AND/OR  DOUBLE 

C  PRECISION  IS  VERY  EASY. 

C 

C  THE  PROGRAM  CALLS  THE  FOLLOWING  SUBROUTINES  WHOSE  SOURCE  CODE  IS 

C  NOT  INCLUDED  HERE: 

C         FFT2C        IMSL  ROUTINE  TO  COMPUTE  FFT  OF  A  COMPLEX 

C  VALUED  SEQUENCE  OF  LENGTH  EQUAL  TO  POWER 

C  OF  TWO 

C  FFT3D         IMSL  ROUTINE  TO  COMPUTE  THE  FFT  OF  A  COMPLEX 

C  VALUED  1,2,  OR  3  DIMENSIONAL  ARRAY 

C         FFTCC        IMSL  ROUTINE  CALLED  BY  FFT3D 

C         MMBSJN       IMSL  ROUTINE  TO  CALCULATE  BESSEL  FUNCTION  OF 

C  THE  FIRST  KIND  OF  ZERO  ORDER  WITH  REAL  ARGUMENT 

C  MMBSIN       IMSL  ROUTINE  TO  CALCULATE  MODIFIED  BESSEL 

C  FUNCTION  OF  THE  FIRST  KIND  OF  ZERO  ORDER  WITH 

C  REAL  ARGUMENT 

C 

C  T.D.M.   20  JANUARY  1987 

Q  it  it  it  it  it  it  it  A  A  A  A  A  tV  A  A  A  A  A  A  A  A  AA  AA  A  A  A  A  A  A  A  iV  A  AAA  it  it  AAA  A  A  A  A  A  A  A  A  *  .V  Vr  rt  A  *  fr  Vr  *■  Yr  iV  iV  rt  Vr  Vr  rt 


C 


REAL  A , A2 , B , B 1 , BASE , BI , BR , C , C2 , CENTER 
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REAL 
REAL 
REAL 
REAL 


HWIDTH , INC , MAXTIME , MEANX , PI 

R123 , RESULT , RH01 , RH02 , RWK , SIGMA , SIGMA2 

TEMPI , TEMP2 , TEMP3 , TEMP4 , TEMP5 

TIME , TS , TWOPI , X , Y , Z , Z2 , ZPRIME , TVECTR 


COMPLEX  C123,CWK,CTEMP1,CTEMP2,CTEMP3 

COMPLEX  PULSE, WORK 

INTEGER  COL.FTYPE.H, I, IER, I JOB, IWK, IWK1, J 

INTEGER  M,N,NAMLEN,NCOLS,NROWS,Rl,R2 

INTEGER  ROW , RHO , SHAPE , STEP , WIDTH , XPOS , YPOS 

DIMENSION  CWK( 64 ) , IWK( 534 ) , IWK1 ( 7 ) ,RH01 ( 64 )  , RWK( 534 ) 

DIMENSION  PULSE(64,64),RESULT(64,64),RH02(64,64) 

DIMENSION  TIME(64),WORK(64,64),ZPRIME(64) 


CHARACTER      FN*8,FNAME*8,T1 ,T2 

CHARACTER      EXT1*4 ,EXT2*4 ,EXT3*4 ,EXT4*4 

CHARACTER      FNAME1*12,  FNAME2*12,  FNAME3*12,  FNAME4*12 

***************************************************************** 

CONSTANTS 
***************************************************************** 

C  SPEED  OF  SOUND  IN  VACUUM  (1500  METERS /SEC) 

C2  SPEED  OF  SOUND  SQUARED  (METERS**2/SEC**2) 

H  HALF  OF  NROWS  (OR  NCOLS)  PLUS  1 

M  INTEGER  USED  TO  SPECIFY  SIZE  OF  VECTOR  SENT  TO 
FFT2C  (REFER  TO  IMSL  SOURCE  CODE  FOR  FFT2C) 

MEANX  USED  FOR  GAUSSIAN  IMPULSE  FUNCTION 

N  USED  ON  BESSEL  FUNCTION  CALLS  TO  SIGNIFY  ZERO  ORDER 

NCOLS  DEFAULT  NUMBER  OF  COLUMNS 

NROWS  DEFAULT  NUMBER  OF  ROWS 

SIGMA  USED  FOR  GAUSSIAN  IMPULSE  FUNCTION 

SIGMA2  SIGMA  SQUARED 

STEP  NUMBER  OF  ROWS  SET  TO  ZERO  -  SIMULATES  STEP  FUNCTION 


TWOPI    2  *  PI 
************************* 

C  =  1500.0 

N  =  1 

M  =  6 

STEP  =  4 

NROWS  =64 

NCOLS  =64 

PI  =  3.14159265 

SIGMA  =1.0 

MEANX  =0.0 

C2  =  C  *  C 

H  =  NROWS/2  +  1 

TWOPI  =  2.0  *  PI 

SIGMA2  =  SIGMA  *  SIGMA 


************************************** 


54 


Q         ***************************************************************** 

C  GET  USER  INPUT 

C     ***************************************************************** 

XPOS  =  1 

YPOS  =  2 

CALL  CLRSCR 

CALL  GOTOXYU.l) 

WRITE(*,*)  'ENTER  REQUESTED  VALUES.   RECOMMENDED  VALUES  IN  (  ) ' 

CALL  GOTOXYC YPOS, XPOS) 

WRITEC*,*)  'ENTER  FILENAME  (8  CHAR  MAX):   * 

CALL  GOTOXYC YPOS, XPOS+3 2) 

READ(*,'(A)')  FN 

CALL  GOTOXY(YPOS+l,XPOS) 

WRITEC*,*)  'SELECT  FILTER  TYPE' 

CALL  GOTOXYC YPOS+2,XPOS+8) 

WRITEC*,*)  '<1>  GP1  CLOSSLESS  MEDIA)' 

CALL  GOTOXYC YPOS+3 , XPOS+8 ) 

WRITEC*,*)  '<2>     GP2  CLOSSY  MEDIA)' 

CALL  GOTOXYCYPOS+4.XPOS) 

WRITEC*, 9)  '  ENTER  1  OR  2:   ' 

READC*,*)  FTYPE 

CALL  GOTOXY(YPOS+5,XPOS) 

WRITEC*,*)  'SELECT  SOURCE  SHAPE' 
'CALL  GOTOXYC YPOS+6 , XPOS+4 ) 

WRITEC*,*)  '<1>   SQUARE' 

CALL  GOTOXYC YPOS+7 .XPOS+4) 

WRITEC*,*)  '<2>  CIRCLE' 

CALL  GOTOXYCYPOS+8, XPOS+4) 

WRITEC*,*)  '<3>  GAUSSIAN' 

CALL  GOTOXYCYPOS+9, XPOS+4) 

WRITEC*,*)  '<4>  PYRAMID' 

CALL  GOTOXYC YPOS+ 10, XPOS+4) 

WRITEC*,*)  '<5>  TRUNCATED  PYRAMID' 

CALL  GOTOXYC YPOS+ 11, XPOS+4) 

WRITEC*,*)  '<6>  5  PULSE  LINEAR  ARRAY  CSQ.)' 

CALL  GOTOXYC YPOS+ 12, XPOS+4) 

WRITEC*,*)  '<7>   5  PULSE  LINEAR  ARRAY  (GAUSS.)' 

CALL  GOTOXYC YPOS+6,XPOS+40) 

WRITEC*,*)  '  <8>   5  PULSE  LINEAR  ARRAY  (STEPPED)' 

CALL  GOTOXYC YPOS+7, XPOS+4 0) 

WRITEC*,*)  '  <9>   5  PULSE  LINEAR  ARRAY  (PARABOLA)' 

CALL  GOTOXY(YPOS+8, XPOS+4 0) 

WRITE (*,*)  '<10>  RESERVED  FOR  FUTURE  USE' 

CALL  GOTOXYC YPOS+9,XPOS+40) 

WRITEC*,*)  '<11>  RESERVED  FOR  FUTURE  USE' 

CALL  GOTOXY(YPOS+10,XPOS+40) 

WRITE (*,*)  '<12>  RESERVED  FOR  FUTURE  USE' 

CALL  GOTOXY(YPOS+ 11, XPOS+4 0) 

WRITE (*,*)  '<13>  RESERVED  FOR  FUTURE  USE' 

CALL  GOTOXY(YPOS+ 12, XPOS+4 0) 

WRITE (*,*)  '<14>  RESERVED  FOR  FUTURE  USE' 
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CALL  G0T0XY(YP0S+13,XP0S) 
WRITE(*,9)  '  ENTER  SELECTION:   * 
READ(*,*)  SHAPE 
IF  (SHAPE  .LT.  6)  THEN 

CALL  GOTOXY(YPOS+14,XPOS) 

WRITE(*,9)  '  ENTER  SOURCE  SIZE  (ODD  INTEGER  <=  63)  (11):   ' 

READ(*,*)  WIDTH 
ENDIF 

CALL  GOTOXY(YPOS+15,XPOS) 
WRITE(*,*)  'ENTER  RHO   (200):  ' 
CALL  GOTOXY(YPOS+15,XPOS+26) 
READ(*.*)  RHO 
CALL  GOTOXY(YPOS+16,XPOS) 
WRITE(*,*)  'ENTER  Z  (METERS)   (0.1):  ' 
CALL  GOTOXY(YPOS+16,XPOS+26) 
READ(*,*)  Z 

CALL  GOTOXY(YPOS+17,XPOS) 
WRITE(*,7)  Z/C 
CALL  GOTOXY(YPOS+17,XPOS+44) 
READ(*,*)  MAXTIME 
IF  (FTYPE  .EQ.  2)  THEN 

CALL  GOTOXY(YPOS+18,XPOS) 

WRITE(*,*)  'ENTER  VALUE  FOR  A   (0.008):   ' 

CALL  GOTOXY(YPOS+18,XPOS+26) 

READ(*,*)  fc 

A2  =  A  *  A 
ENDIF 
C 

7     FORMATC  ENTER  MAXTIME  (REAL  >  '  ,  1PEU  .  6,  '  ) :  '  ) 
9     FORMAT (A\) 

HWIDTH  =  REAL(WIDTH)/2 
Z2  =  Z  *  Z 
TS  =  Z/C 
CALL  CLRSCR 
CALL  GOTOXY(l,l) 
WRITE(*,*)  'CALCULATING...' 
C 

C  INITIALIZE  ARRAYS  AND  VECTORS  TO  ZERO 

Q  ***************************************************************** 

DO  5  ROW  =  l.NROWS 
TIME(ROW)  =0.0 
RHOKROW)  =0.0 
DO  5  COL  =  l.NCOLS 

PULSE (COL, ROW)  =  (0.0,0.0) 
WORK(COL.ROW)  =  (0.0,0.0) 
RESULT (COL, ROW)  =0.0 
RH02(COL,ROW)  =  0.0 
5     CONTINUE 
C 
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c 


c  ***************************************************************** 

C  OPEN  OUTPUT  FILES 

***************************************************************** 
C        PROGRAM  GETS  NAME  OF  PROBLEM  <FN>  FROM  USER  AND  FORMS: 
C           FILENAME. EXT           FILE  DESCRIPTION 
c  

C  <FN>.IMP     ARRAY  OF  INPUT  IMPULSE  FUNCTION 

C  <FN>.TXT     SUMMARY  INFORMATION  FOR  THIS  PROBLEM 

C  <FN>.DAT     FINAL  RESULTS  ARRAY 

C 

C        NOTE  1:   ALL  FILES  ARE  IN  ASCII  TEXT  FORMAT 

n  ***************************************************************** 

EXT1  =  ' .IMP' 
EXT2  =  ' . TXT ' 
EXT3  =  ' . DAT ' 
NAMLEN  =  0 
Tl  =  '  ' 
DO  10  1=1,8 

T2  =  FN(I:I) 

IF  (T2  .NE.  Tl)  THEN 
FNAME(I:I)  =  FN(I:I) 
NAMLEN  =  NAMLEN  +  1 

ENDIF 
10    CONTINUE 

FNAMEK1:  NAMLEN)  =  FNAME 
FNAMEl(NAMLEN+l:NAMLEN+4)  =  EXT1 
FNAME2C1: NAMLEN)  =  FNAME 
FNAME2(NAMLEN+1:NAMLEN+A)  =  EXT2 
FNAME3  ( 1 :  NAMLEN )  =  FNAME 
FNAME3(NAMLEN+l:NAMLEN+4)  =  EXT3 
IF  (NAMLEN  .LT.  8)  THEN 

FNAME1(NAMLEN+5:12)  =  ' 

FNAME2(NAMLEN+5:12)  =  '         ' 

FNAME3(NAMLEN+5:12)  =  ' 
ENDIF 

OPEN (1,FILE=FNAME1.STATUS=' NEW* ) 
OPEN ( 2 , F ILE=FNAME2 , STATUS= ' NEW  * ) 
OPEN(3,FILE=FNAME3,STATUS='NEW' ) 
C 

Q        ***************************************************************** 

C  WRITE  HEADER  INFORMATION  TO  <FN>.TXT 

Q  ***************************************************************** 

WRITE(2,11)  FN 

IF  (FTYPE  .EQ.  1)  THEN 

WRITEC2,*)  '   FILTER  IS  GP1  FOR  LOSSLESS  MEDIA' 
ELSEIF  (FTYPE  .EQ.  2)  THEN 

WRITE(2,*)  '   FILTER  IS  GP2  FOR  LOSSY  (LINEAR)  MEDIA* 
ENDIF 
IF  (SHAPE  .EQ.  1)  THEN 

WRITE (2, 12)  WIDTH, NCOLS.NROWS 


57 


ELSEIF  (SHAPE  .EQ.  2)  THEN 

WRITE (2, 13)  WIDTH, NCOLS.NROWS 
ELSEIF  (SHAPE  .EQ.  3)  THEN 

WRITE (2, 14)  WIDTH, NCOLS.NROWS 
ELSEIF  (SHAPE  .EQ.  4)  THEN 

WRITE (2, 111)  WIDTH , NCOLS , NROWS 
ELSEIF  (SHAPE  .EQ.  5)  THEN 

WRITE (2, 112)  WIDTH, NCOLS, NROWS 
ELSEIF  (SHAPE  .EQ.  6)  THEN 

WRITE(2,113) 
ELSEIF  (SHAPE  .EQ.  7)  THEN 

WRITE (2, 114) 
ELSEIF  (SHAPE  .EQ.  8)  THEN 

WRITE(2,115) 
ELSEIF  (SHAPE  .EQ.  9)  THEN 

WRITE(2,116) 
ENDIF 

WRITE(2,15)  STEP 
WRITE(2,16)  Z 
WRITE(2,17)  RHO 
WRITE(2,18)  MAXTIME 
IF  (FTYPE  .EQ.  2)  THEN 

WRITE(2,19)  A 
ENDIF 

SUMMARY  INFORMATION  FOR  \A/) 

IMPULSE  FUNCTION  IS  A  SQUARE  WITH  WIDTH  ',12,'  ON  A  BAS 
'  X* ,13) 

IMPULSE  FUNCTION  IS  A  CIRCLE  WITH  DIAMETER  ',12,'  ON  A 
,13, '  X' ,13) 

IMPULSE  FUNCTION  IS  GAUSSIAN  WITH  DIAMETER  ',12,*  ON  A 
,13,'  X' ,13) 

IMPULSE  FUNCTION  IS  PYRAMID  WITH  WIDTH  ',12,'  ON  A  BAS 
,*  X' ,13) 

IMPULSE  FUNCTION  IS  TRUNCATED  PYRAMID  WITH  WIDTH  ' , 12 , ' 
■+ON  A  BASE  OF' ,13, '  X' ,13) 
113    FORMATC    IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  SQUARE  PU 


IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  GAUSSIAN 
IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  STEPPED  A 
IMPULSE  FUNCTION  IS  5  PULSE  LINEAR  ARRAY  WITH  PARABOLIC 


11 

FORMATC : 

12 

FORMATC 

• 

+E  OF*, 13 

13 

FORMAT ( ' 

+BASE  OF' 

14 

FORMATC 

+BASE  OF' 

111 

FORMATC 

+E  OF' ,13 

112 

FORMATC 

+LSES '  ) 

114 

FORMAT ( ' 

+PULSES '  ) 

115 

FORMATC 

+MPLITUDE* ) 

116 

FORMATC 

+  AMPLITUDE 

15 

FORMATC 

16 

FORMATC 

17 

FORMATC 

18 

FORMATC 

19 

FORMATC 

C 

STEPSIZE  =  ' ,12) 

Z  =  ' .1PE14.7,'  METERS') 

RHO  =  ' ,13) 

MAXTIME  =  \1PE14.7,'  SECONDS') 

A  =  ' .1PE14.7) 
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c  ***************************************************************** 

C  INITIALIZE  TIME  VECTOR 

C     ***************************************************************** 

INC  =  (MAXTIME-TS)/REAL(NROWS-STEP-l) 
DO  20  I  =  (STEP+l),NROWS 

TIME(I)  =  TS  +  (I-STEP-1)  *  INC 

20  CONTINUE 
C 

Q  ***************************************************************** 

C  OUTPUT  TIME  SUMMARY  TO  <FN>.TXT 

Q        ***************************************************************** 

WRITE(2,*) 

WRITE(2,9)  'TIME  SUMMARY:' 
WRITE(2,21)  TS.STEP+1 
WRITE (2 ,22)  INC 

21  FORMATC    ZERO  TIME  STARTS  AT  T=Z/C  =',1PE14.7,'  SECONDS  IN  ROW( ' 
+,12,')') 

22  FORMATC      AND  INCREMENTS  BY  '.1PEU.7,'  SECONDS  PER  ROW'/) 
C 

Q  ***************************************************************** 

C         CREATE  Z- PRIME  VECTOR:   SQUARE  EACH  VALUE  OF  TIME, 
C  MULTIPLY  BY  C  SQUARED  AND  SUBTRACT  Z  SQUARED.  THEN 

C  TAKE  SQUARE  ROOT. 

Q  ***************************************************************** 

DO  30  I  =  (STEP+2),NROWS 

ZPRIME(I)    =  SQRTUTIME(I)    *   TIME(I)    *   C2)    -    Z2) 
30  CONTINUE 

C 

Q  ************************************************* 

C  INITIALIZE  RHO  VECTOR  AND  GENERATE  33  X  64  RHO  MATRIX 

C       •   ONLY  ONE  HALF  OF  THE  MATRIX  HAS  TO  BE  CALCULATED  SINCE 
C  THE  OTHER  SIDE  IS  SYMMETRICAL 

Q  ***************************************************************** 

INC  =■  REAL(RHO)/REAL(NCOLS/2-l) 
BASE  =  REAL(-RHO)  -  INC 
'RHOl(l)  =  BASE 
DO  50  I  =  l.NCOLS-1 

RHOKI+1)  =  BASE  +  I  *  INC 
50    CONTINUE 
C 

DO  60  I  =  l.NROWS 

RH02(H,I)  =  RHOl(I) 
RH02(I,H)  =  RHOl(I) 
60    CONTINUE 
C 

DO  70  ROW  =  l.NROWS 
DO  70  COL  =  1,H 

RH02(COL,ROW)  =  SQRT(RH02(COL,H)**2  +  RH02(H,ROW)**2) 
70    CONTINUE 
C 
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Q  ***************************************************************** 

C  OUTPUT  RHO  SUMMARY  TO  <FN>.TXT 

Q  ***************************************************************** 

WRITE(2,9)  'RHO  SUMMARY:' 
WRITE(2,44)  H.NCOLS/2+1 
WRITE(2,A1)  BASE 
WRITE(2.42)  INC 
WRITE(2,43)  RH02(1,1) 

41  FORMATC      BASE  VALUE  =  (-RHO  -  INC)     =  \1PE14.7) 

42  FORMATC      INCREMENT  =  RHO/ (64/2  -  1)    =  ' .1PE14.7) 

43  FORMATC      MAXIMUM  VALUE  IS  AT  RH02(1,1)  =  \1PE14.7) 

44  FORMATC    ARRAY  IS  CENTERED  ABOUT  C , 13 , ' , ',13 , ' ) ' ) 
C 

Q  ***************************************************************** 

C         GENERATE  EXCITATION  -  CENTER  POINT  MUST  BE  AT  (33,33) 

Q  ***************************************************************** 

C 

C     ***   GENERATE  SQUARE  EXCITATION   *** 
IF  (SHAPE  .EQ.  1)  THEN 

DO  1200  ROW  =  (NROWS/2+1-WIDTH/2), (NROWS/2+1+WIDTH/2) 
DO  1200  COL  =  (NC0LS/2+1-WIDTH/2), (NC0LS/2+1+WIDTH/2) 
PULSE (COL, ROW)  =  (1.0,0.0) 
1200     CONTINUE 
C 

C     ***   GENERATE  CIRCULAR  PULSE   *** 
•   ELSEIF  (SHAPE  .EQ.  2)  THEN 

DO  1300  ROW  =  (NROWS/2+l-WIDTH/2),(NROWS/2+l+WIDTH/2) 
DO  1300  COL  =  (NCOLS/2+1-WIDTH/2), (NCOLS/2+1+WIDTH/2) 

IF  (SQRT(REAL(COL-H)**2+REAL(ROW-H)**2)  XT.  HWIDTH)  THEN 

PULSE (COL, ROW)  =  (1.0,0.0) 
END  IF 
1300     CONTINUE 
C 

C     ***   GENERATE  CIRCULAR  EXCITATION  WITH  GAUSSIAN  AMPLITUDE   *** 
ELSEIF  (SHAPE  .EQ.  3)  THEN 

TEMPI  =  1/SQRT(TW0PI*SIGMA2) 
TEMP2  =  1/(2*SIGMA2) 

DO  1400  ROW  =  (NROWS/2+l-WIDTH/2),(NROWS/2+l+WIDTH/2) 
DO  1400  COL  =  (NC0LS/2+1-WIDTH/2), (NCOLS/2+1+WIDTH/2) 
TEMP3  =  SQRT(REAL(COL-H)**2+REAL(ROW-H)**2) 
IF  (TEMP3  .LT.  HWIDTH)  THEN 

PULSE (COL, ROW)  =  TEMP1*EXP(-TEMP2*(TEMP3-MEANX)**2) 
ENDIF 
1400     CONTINUE 
C 


***   GENERATE  PYRAMID  EXCITATION 
ELSEIF  (SHAPE  .EQ.  4)  THEN 
Rl  =  (NROWS/2)  -  (WIDTH/2) 
R2  =  (NROWS/2  +  2)  +  (WIDTH/2) 
INC  =  1.0/ (WIDTH/2  +  1) 
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DO  1500  I  =  1, WIDTH/2 

DO  1500  ROW  =  R1+I.R2-I 
DO  1500  COL  =  R1+I.R2-I 

PULSE (COL, ROW)  =  INC  *  I 
1500     CONTINUE 

PULSE(H.H)  =  1.0 
C 

C     ***   GENERATE  TRUNCATED  PYRAMID  EXCITATION   *** 
ELSEIF  (SHAPE  .EQ.  5)  THEN 
Rl  =  (NROWS/2)  -  (WIDTH/2) 
R2  *   (NROWS/2  +  2)  +  (WIDTH/2) 
INC  =  1.0/(NCOLS/2) 

BASE  =  (NCOLS/2  -  WIDTH/2  +  1)  *  INC 
INC  -  1.0/H 
DO  1600  I  =  1, WIDTH/2 

DO  1600  ROW  =  R1+I.R2-I 
DO  1600  COL  =  R1+I.R2-I 

PULSE (COL, ROW)  =  BASE  +  INC  *  I 
1600     CONTINUE 

PULSE(H.H)  =  1.0 
C 

C     ***   GENERATE  5  ELEMENT  LINEAR  ARRAY  -  SQUARE   *** 
ELSEIF  (SHAPE  .EQ.  6)  THEN  . 
WIDTH  =  5 

DO  1700  I  =  1, WIDTH 
COL  =  I  *  10 
DO  1700  J  =  1, WIDTH 
COL  =  COL  +  1 

DO  1700  ROW  =  (H  -  WIDTH/2), (H  +  WIDTH/2) 
PULSE (COL, ROW)  =  (1.0,0.0) 
1700     CONTINUE 
C 

C     ***   GENERATE  5  ELEMENT  LINEAR  ARRAY  -  GAUSSIAN   *** 
ELSEIF  (SHAPE  .EQ.  7)  THEN 
WIDTH  =  5 

TEMPI  =  1/SQRT(TW0PI*SIGMA2) 
TEMP2  =  1/(2*SIGMA2) 

DO  1800  ROW  =  (NROWS/2+l-WIDTH/2),(NROWS/2+l+WIDTH/2) 
DO  1800  COL  =  (NCOLS/2+l-WIDTH/2),(NCOLS/2+l+WIDTH/2) 
TEMP3  =  SQRT(REAL(COL-H)**2+REAL(ROW-H)**2) 
IF  (TEMP  .LT.  REAL (WIDTH/2))  THEN 

TEMPA  =  TEMP1*EXP(-TEMP2*(TEMP3-MEANX)**2) 
PULSE (COL-20, ROW)  =  TEMPA 
PULSE (COL -10, ROW)  =  TEMPA 
PULSE (COL, ROW)     =  TEMPA 
PULSE(COL+10,ROW)  =  TEMPA 
PULSE (COL+20, ROW)  =  TEMPA 
ENDIF 
1800     CONTINUE 
C 
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C  ***        GEKERATE   5  ELEMENT  LINEAR  ARRAY   -   STEPPED        *** 

ELSEIF    (SHAPE    .EQ.    8)   THEN 
WIDTH  =   5 
TEMPI  =   1.0/3.0 
TEMP2  =   2.0/3.0 
TEMP3   =1.0 
DO   1900   COL   =    1, WIDTH 

DO  1900  ROW  =  (H-WIDTH/2),(H+WIDTH/2) 
PULSE (COL+ 10, ROW)  =  TEMPI 
PULSE(COL+20,ROW)  =  TEMP2 
PULSE (COL+30, ROW)  =  TEMP3 
PULSE(COL+40,ROW)  =  TEMP2 
PULSE(COL+50,ROW)  =  TEMPI 
1900  CONTINUE 

C 

C  ***   LINEAR  ARRAY  OF    5  ELEMENTS   -    PARABOLIC  AMPLITUDE        *** 

C  IMPLEMENTS   THE   FOLLOWING  FORMULA: 

C  (X   -   H)**2   =   -4    *   P   *    (Y   -   K) 

C  WITH 

C  H  =   0 

C  P   =   256 

C  K  =   256 

C  TO  GIVE  A  DOWNWARD  OPENING   PARABOLA  THAT   HAS   A  PEAK  OF    1 . 0 

C  AND    IS   EQUAL   TO   ZERO  AT  X  =   -32,    32    (ROWS   2   AND   64). 

C  EQUATION   IS   EVALUATED   AT  X  =   0 ,    10   AND   20. 

C 

ELSEIF    (SHAPE    .EQ.    9)    THEN 
WIDTH   =■   5 

TEMPI  =  (A.*256.  -  20**2)/  (4  .*256  .  ) 
TEMP2  =  (4.*256.  -  10**2)/ (4 .*256. ) 
TEMP3  =  (4.*256.  -  0**2)/ ( A .*256. ) 
DO   2000   COL   =    1, WIDTH 

DO  2000  ROW  =  (H-WIDTH/2),(H+WIDTH/2) 
PULSE (COL+ 10, ROW)  =  TEMPI 
PULSE(COL+20,ROW)  =  TEMP2 
PULSE (COL+30, ROW)  =  TEMP3 
PULSE(COL+40,ROW)  =  TEMP2 
PULSE(COL+50,ROW)  =  TEMPI 
2000  CONTINUE 

C 

ENDIF 
C 

C  OUTPUT    IMPULSE   FUNCTION  ARRAY  TO  <FN>.IMP 

Q  *tt**»***t***iV*******»*»****4*»*«*(r**ilil*M»*»«*t*4**H4t**»****** 

DO   3000    I   =    l.NROWS 

WRITE(1,102)    (REAL(PULSE(J,I)),J  =    l.NCOLS) 
3000      CONTINUE 

CLOSE ( 1 ) 
C 
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c  ***************************************************************** 

C        PASS  PULSE  TO  2-D  FFT  SUBROUTINE.   I JOB  =  1  MEANS  TAKE  FFT 

C  ***************************************************************** 

IJOB  =  1 

CALL  FFT3D ( PULSE , NCOLS , NROWS , NCOLS , NROWS , 1 , I JOB , IWK , RWK , CWK ) 
C 

Q  ***************************************************************** 

C       SHIFT  RESULTS  -  IE:  SWITCH  QUADRANTS  1  AND  3,  2  AND  A  IN  THE 

C  COL, ROW  PLANE 

r     ***************************************************************** 

CALL  SHIFT  ( PULSE , NROWS , NCOLS ) 
C 

Q  ***************************************************************** 

C  MAIN  PROGRAM  STARTS  HERE 

C        CALCULATIONS  START  AT  ROW  =  STEP+1. 

C        NOTE  THAT  ONLY  ONE  QUARTER  OF  THE  WORK  MATRIX  HAS  TO 

C        BE  CALCULATED.   THE  REMAINING  THREE  QUARTERS  IS  FORMED  BY 

C        FOLDING  AND  FLIPPING.   THIS  IS  POSSIBLE  DUE  TO  SYMMETRY. 

C  ***************************************************************** 

CALL  GOTOXY(l.l) 

WRITE (*,*)  'CALCULATING  ROW   ' 

IF  (FTYPE  .EQ.  1)  THEN 

DO  200  I  =  STEP+1, NROWS 
CALL  GOTOXY(l,18) 
WRITE(*,101)  I 

DO  210  ROW  =  1,H 
DO  210  COL  =  1,H 

Bl  =  ZPRIME(I)  *  RH02(COL,ROW) 
CALL  MMBS JN(B1,N, TEMPI, IER) 
WORK(COL,ROW)  =  TEMPI  *  PULSE ( COL , ROW ) 
210         CONTINUE 
C 
C     ***************************************************************** 

C  FLIP  QUADRANT  2  TO  QUADRANT  3 

Q       ***************************************************************** 

J  =  0 

DO  220  ROW  =  H+l,  NROWS 
J  =  J  +  2 
DO  220  COL  =  1,H 

WORKCCOL.ROW)  =  WORK (COL, ROW- J) 
220         CONTINUE 
C 

Q  ***************************************************************** 

C  FLIP  LEFT  SIDE  TO  RIGHT  SIDE 

Q  ***************************************************************** 

J  =  0 

DO  215  COL  =  H+l, NCOLS 
J  =  J  +  2 
DO  215  ROW  =  1, NROWS 

WORK(COL,ROW)  =  WORK (COL -J, ROW) 
215         CONTINUE 


63 


Q  ***************************************************************** 

C  TAKE  INVERSE  FFT  OF  EACH  COLUMN.   THEN  TAKE  INVERSE  FFT  OF 

C  ROW  33  ONLY,  SINCE  THIS  IS  THE  ROW  CONTAINING  CENTRAL 

C  INFORMATION.   ROW  33  BECOMES  THE  NEXT  ROW  OF  THE  FINAL 

C  OUTPUT  MATRIX. 

Q  ***************************************************************** 

DO  A 15  COL  =  1, NCOLS 
DO  405  ROW  =  1, NROWS 

CWK(ROW)  =  CON JG( WORK (COL, ROW)) 
405  CONTINUE 

CALL  FFT2C  (CWK,M,IWK1) 
DO  410  ROW  =  l.NROWS 

WORK(COL.ROW)  =  CWK(ROW) 
410  CONTINUE 

415         CONTINUE 
C 

DO  440  COL  =  1, NCOLS 

CWK(COL)  =  WORK(COL.H) 
440         CONTINUE 

CALL  FFT2C  (CWK.M.IWK1) 
C 

R123  =  NROWS  *  NCOLS 
C123  =  CMPLX(R123,0.0) 
DO  450  COL  =  1, NCOLS 

RESULT(COL.I)  =  ABS(REAL( (CONJG(CWK(COL) ) /C123) ) ) 
450         CONTINUE 
C 

200      CONTINUE 
C 

ELSEIF(FTYPE  .EQ.  2)  THEN 
TEMPI  =  A  *  C  /  2 
TEMP2  =  A2  *  C2  /  4 
CALL  GOTOXY(l,l) 
WRITE (*,*)  'CALCULATING  ROW   ' 
DO  300  I  =  STEP+1, NROWS 
CALL  GOTOXY(l,18) 
WRITE(*,101)  I 

DO  310  ROW  =  1,H 
DO  310  COL  =  1,H 

TEMP3  =  ZPRIME(I)  *  SQRT(ABS(RH02(COL,ROW)**2-TEMP2) ) 
IF  (ABS(RH02(COL,ROW))  .GT.  TEMPI)  THEN 

CALL  MMBSJN(TEMP3,N,TEMP4,IER) 
ELSE 

CALL  MMBSIN(TEMP3,N,TEMP4,IER) 
END  IF 

WORK(COL.ROW)  =  TEMP4  *  EXP(-A*C2*TIME(I ) )  * 
&  PULSE (COL, ROW) 

310         CONTINUE 
C 


64 


e    ***************************************************************** 

C  FLIP  QUADRANT  2  TO  QUADRANT  3 

Q  ***************************************************************** 

J  =  0 

DO  320  ROW  =  H+l.NROWS 
J  =  J  +  2 
DO  320  COL  =  l.H 

WORK(COL,ROW)  =  WORK (COL, ROW- J) 
320         CONTINUE 

C 

C     ***************************************************************** 

C  FLIP  LEFT  SIDE  TO  RIGHT  SIDE 

Q         ***************************************************************** 

J  =  0 

DO  315  COL  =  H+l, NCOLS 
J  =  J  +  2 
DO  315  ROW  =  l.NROWS 

WORK(COL,ROW)  =  WORK (COL -J, ROW) 
315         CONTINUE 
C 
C     ***************************************************************** 

C  TAKE  INVERSE  FFT  OF  EACH  COLUMN.   THEN  TAKE  INVERSE  FFT  OF 

C  ROW  33  ONLY,  SINCE  THIS  IS  THE  ROW  CONTAINING  CENTRAL 

C  INFORMATION.   ROW  33  BECOMES  THE  NEXT  ROW  OF  THE  FINAL 

C  OUTPUT  MATRIX. 

Q  ***************************************************************** 

DO  615  COL  =  1, NCOLS 
DO  605  ROW  =  l.NROWS 

CWK(ROW)  =  CON JG( WORK (COL, ROW)) 
605  CONTINUE 

CALL  FFT2C  (CWK.M.IWK1) 
DO  610  ROW  =  l.NROWS 

WORK(COL,ROW)  =  CWK(ROW) 
610  CONTINUE 

615         CONTINUE 
C 

DO  640  COL  =  l.NCOLS 

CWK(COL)  =  WORK(COL.H) 
640         CONTINUE 

CALL  FFT2C  (CWK,M,IWK1) 
C 

R123  =  NROWS  *  NCOLS 
C123  =  CMPLX(R123,0.0) 
DO  650  COL  =  1, NCOLS 

RESULT(COL.I)  =  ABS(REAL( (CONJG(CWK(COL) ) /C123) ) ) 
650         CONTINUE 
C 

300      CONTINUE 
C 

ENDIF 
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500 

C 

C 

C 

C 


510 
C 

101 
102 

c 
c 
c 
c 
c 
c 


***************************************************************** 

SET  ROWS  1  TO  STEP  TO  ZERO  TO  SIMULATE  STEP  FUNCTION 
***************************************************************** 

DO  500  ROW  =  l.STEP 

DO  500  COL  =  l.NCOLS 

RESULT (COL, ROW)  =0.0 
CONTINUE 

***************************************************************** 

OUTPUT  FINAL  ARRAY  TO  DISK 
***************************************************************** 

DO  510  I  =  l.NROWS 

WRITEO.102)    (RESULT(JrI),J  =   l.NCOLS) 
CONTINUE 

FORMAT  (12) 
FORMAT  (64F16.7) 


*** ************************************************************** 

CLOSE  REMAINING  OPEN  FILES 
***************************************************************** 

CLOSE (2) 
CLOSEO) 

***************************************************************** 

REMIND  USER  OF  DATA  FILES  CREATED 
***************************************************************** 


CALL  CLRSCR 
WRITE(*,*)  'FINISHED.* 


WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

*) 

WRITE (* 

") 

WRITEC* 

*) 

WRITE (* 

*) 

STOP 

END 

*)  '    THE  FOLLOWING  ASCII  TEXT  FILES  ARE  IN  THE  DEFAULT' 
DIRECTORY:' 


FILENAME 


DESCRIPTION' 


'.FNAMEl,'    INPUT  IMPULSE  FUNCTION' 
',FNAME2,'    SUMMARY  INFORMATION' 
' , FNAME3 , '    OUTPUT  ARRAY ' 
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c  ****************************************************************** 

c 

C  *  SUBROUTINE  GOTOXY  * 

C  *     MOVES  CURSOR  TO  LINE  X,  COLUMN  Y  * 

C  *     CALL:   CALL  GOTOXY(X.Y)  * 


****************************************************************** 

SUBROUTINE  GOTOXY (X,Y) 
CHARACTER*1  CI ,C2,C5 ,C8,LC(5) 
CHARACTER* 5  CBUFF 
INTEGER*2  IC(4),L,X,Y 

EQUIVALENCE  (CI, IC(1) ) , (C2, IC(2) ) , (C5, IC(3) ) , (C8, IC( A) ) , 
+  (CBUFF,LC(D) 
DATA  IC/16#1B,16#5B,16#3B, 16*66/ 


L  =  10000+100*X+Y 
C 
C  ***  WRITE  ESCAPE  CODES  TO  CBARACTER  BUFFER 

WRITE (CBUFF, 2)  L 
2     FORMAT(I5) 
C 
C  ***  WRITE  ESCAPE  CODES  TO  TBE  DISPLAY 

WRITE (*,1)  C1,C2,LC(2),LC(3),C5,LC(4),LC(5),C8 
1    CF0RMAT(1X,8A1,\) 
C 
C  ***  END  OF  SUBROUTINE 

RETURN 

END 
C 

Q  ***************************************************************** 

C  *  * 

C  *  SUBROUTINE  CLRSCR                                       * 

C  *  SUBROUTINE  TO  CLEAR  THE  DISPLAY                        * 

C  *  CALL:   CALL  CLRSCR                                   * 

C  *  * 

Q  ***************************************************************** 

c 

SUBROUTINE  CLRSCR 

CHARACTERS  C1,C2,C3,C4 

INTEGER*2  IC(A) 

EQUIVALENCE  (CI,  IC(1) ) ,  (C2,  IC(2) ) ,  (C3  ,  IC(3) ) ,  (C4  ,  IC( <* ) ) 

DATA  IC/ 16#1B , 16#5B , 16#32 , 16#4A/ 
C 
C  ***  WRITE  ESCAPE  CODE  TO  DISPLAY 

WRITE(*,1)  C1,C2,C3,C4 
1     FORMAT ( IX, 4A1) 
C 
C  ***  END  OF  SUBROUTINE 

RETURN 

END 
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c 

Q  ***************************************************************** 

C  *  SUBROUTINE  TO  SHIFT  A  SQUARE  MATRIX  -  EXCHANGES  QUADRANTS  ONE  * 

C  *  AND  THREE,  AND  TWO  AND  FOUR.  * 

C  *  USE:   CALL  SHIFT (X.NROWS.NCOLS)  * 

C  *       WHERE  X  IS  A  SQUARE  ARRAY  -  REAL  OR  COMPLEX  * 

C  *       NROWS.NCOLS  IS  THE  ARRAY  SIZE  * 

C  *  * 

Q  ***************************************************************** 

SUBROUTINE  SHIFT  (X.NROWS.NCOLS) 
C 
C 

INTEGER  NROWS.NCOLS, ROW, COL ,R2,C2 
COMPLEX  X(NC0LS,NR0WS),T1,T2 
C 

R2  =  NROWS/2 
C2  =  NCOLS/2 
C 

DO  10  ROW  =  1,R2 
DO  10  COL  =  1.C2 
Tl  =  X(COL.ROW) 
X(COL.ROW)  =  X(COL+C2,ROW+R2) 
X(C0L+C2,R0W+R2)  =  Tl 
T2  -  X(C0L+C2,R0W) 
X(COL+C2,R'>J)  =  X(C0L,R0W+R2) 
X(C0L,R0W+R2)  =  T2 
10    CONTINUE 
RETURN 
END 
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